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We model an interface layer connecting two parts of a solid body by N parallel elastic springs 
connecting two rigid blocks. We load the system by a shear force acting on the top side. The 
springs have equal stiffness but are ruptured randomly when the load reaches a critical value. For 
^vq \ the considered system, we calculate the shear modulus, G, as a function of the order parameter, <fi, 

describing the state of damage, and also the "spalled" material (burst) size distribution. In particular, 
we evaluate the relation between the damage parameter and the applied force and explore the 
behaviour in the vicinity of material breakdown. Using this simple model for material breakdown, 
we show that damage, caused by applied shear forces, is analogous to a first-order phase transition. 
The scaling behaviour of G with <f> is explored analytically and numerically, close to 4> = and <f> — 1 
ryj ' and in the vicinity of <j> c , when the shear load is close but below the threshold force that causes 

material breakdown. Our model calculation represents a first approximation of a system subject to 
wear induced loads. 
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The damage of solids subject to external loads is a progressive physical process, which eventually may lead to mate- 
rial breakdown. At the microscopic scale, damage is caused by the localisation of the stress field in the neighbourhood 
of defects or at interfaces that induce the breaking of the atomic bonds. At a mesoscopic scale there is the interaction 
and coalescence of microscopic cracks or vacancies (voids) , which together initiate a macroscopic crack. 

One category of material damage problems, of great technological and fundamental interest, is the deterioration of 
solid surfaces in contact induced by relative motion. These problems belong to the field of tribology, in which friction 
and wear of materials are studied Q, Q • The first phenomenological approach to the problem of wear dates back to 
the 1940's, when Holm Q provided a simple description for wear of solids, which later was re-derived and applied 
by Archard [1] . The Holm- Archard wear description states that the amount of wear or the worn mass between two 
rubbing surfaces is proportional to sliding work (or dissipation energy due to friction) divided by the hardness of 
the surface; simply m oc F^x/p, where m is the worn mass, Fn the normal force or the force applied normal to the 
interface, x the distance slid and p is the penetration hardness, which is related to the surface tension of the solid, 
7, in the manner p oc 7 3 , see 0. A mass loss from a solid body associates with heat dissipated to the environment, 
or from thermodynamics, dQ < TdS, where dQ is the increment of heat exchange, T the temperature and dS is the 
change in entropy. A shear force, Fs, applied to one of the bodies in a direction parallel to the interface plane, causes 
slip (a relative motion between parts of the interface) when Fs = /i-F/v, where /i is the coefficient of friction. This is 
Amontons's second law of macroscopic friction, formulated more than 200 years ago Q. Setting dQ — [iFiydx and 
fi — G/p, with G denoting the shear strength of softer material, we can express the Holm-Archard description as 
q , dm < TdS/G, i.e. the worn mass is proportional to the production of entropy and inversely to the shear strength of 
■ softer material. 

The aim of our study is to investigate the manner in which a solid material breaks down, when it is subjected to 
a shear force that reaches a certain critical value, by using a simple model for material damage. The system under 
study is described by N parallel elastic springs connecting two rigid blocks. The system is loaded with a shear force 
acting on the top side. If the relative motion is slow and the system is well cooled, it is possible to neglect the 
heat generated at the contact. Furthermore, neglecting the influence of the normal load, which mainly introduces 
compressive stresses, and supposing that the crack path is limited to the interface between the solid blocks, we arrive 
at the system studied in this paper. The model can represent rupture of a weld between two solid bodies due to shear 
forces. 

A solid body subject to stress can be considered as being in a metastable state. It can transform in a self-organised 
manner to stable fracture state by formation of cracks. This self-organised behaviour of de-bonding in solids, which is 
the start of the damage process, is theoretically analogous the phenomenon of nucleation in first-order phase transition 
HHQ. Zapperi et al. @, @ have utilised this analogy to study the global breakdown of disordered media under 
external loads. They used a two-dimensional lattice model, in which each bond of the lattice represents an elastic 
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spring that breaks when it is stretched beyond a threshold value governed by a probability distribution. They studied 
numerically the random fuse model (RFM) and a spring network in the framework of a mean-field theory. They 
observed that the breakdown is preceded by avalanche events, which is analogous to the formation of the droplets 
observed near a spinodal decomposition flfl ITll | . 

Lattice or discrete modelling approach to material damage dates back to the 1926 work of Peirce, what is presently 
known as the global sharing fibre bundle 0|. Peirce considered N parallel fibres, each with its own rupture threshold. 
The fibres were linked in such a way that when one link failed, the load was transformed to and shared equally among 
all the intact links. Later, Daniels studied and extended this model in detail and determined the asymptotic 
distribution of the bundle strength for large N. Sornette evaluated the failure characteristic of N independent 
vertical lines linked in parallel with identical spring constants and random failure thresholds. He showed that the 
rupture properties of this system is quite distinct from systems that are linked in series. Furthermore, Sornette 
identified different regimes of rupture depending on the range of stresses, (i) An elastic-reversible behaviour was 
observed for stresses below certain threshold, o\ at which the first link failed, (ii) a stress range a\ < o~n , where the 
system exhibited a non-linear elastic but reversible behaviour, and finally, (iii) for a > <tn global failure occurred. 
Moreover, Sornette showed that the asymptotic properties of global failure threshold of the network system can be 
explained in terms of the central limit theorem of statistics. A review of rupture models can be found in Sornette's book 
The book by Hermann and Roux offers a review of the field up to 1990 ^f|. Hemmer and Hansen studied and 
calculated the distribution of burst avalanches in fibre bundles. They found that for a large class of failure threshold 
distributions, the bursts are distributed according to an asymptotic scaling law s -5 / 3 for simultaneous rupture of s 
fibres. 

In this paper, we pursue a similar approach as Zapperi et al. ||| to investigate analytically and numerically the 
breakdown of disordered media subject to surface loads. We envisage that the solid, subjected to shear loading, 
contains an interface layer, which consists of a network of parallel springs. We study the behaviour of the shear 
modulus, the order parameter characterising the state of damage and burst (avalanche) size distribution. We compare 
the universality of our results with some other failure models, e.g. the fibre bundle models and discuss our results in 
the context of the magnetic system. 

The organisation of this paper is as follows: In section |Hj we present the physical model and its associating con- 
stitutive relations. In section Hill we outline the method of computation, numerical and analytical, for the evaluation 
of the shear modulus, the variation of the applied load with the order parameter, the critical exponents and the 
avalanche size distribution. The results of our computations are presented in section IWl and discussed in section Fvl 
Some relating background material is placed in the appendix. 



II. MODEL 



Consider an interface layer connecting two parts of a solid body consisting of N equi-spaced parallel springs of equal 
length, I, with spring constants, kj, and randomly distributed rupture deformation thresholds, Dj, where j £ [1,N]. 
Initially the springs constants will be set equal (kj — k). If a spring is stretched beyond its rupture threshold, it will 
loose its load carrying capacity, i.e., kj — 0. The springs are loaded with a force, F, acting parallel to the upper 
surface simulating the shear component of a wear-inducing load, Fig. Each spring carries the load, Fj — kjAj, 
where Aj =lj — I denotes the displacement of spring j. Assuming linear deformation, the deformation of spring i can 
be written as 

N -i A i- 1 A 
A, = — A 1 + — A„, (1) 

where Ai and A at is the displacement of the end springs. The state of equilibrium gives 



N N N 

i=l i=l i=2 

where the equation on the right states that the sum of the moments of all the forces acting on the body is zero and 
a denotes the lattice spacing. Equations and l|2j are combined to find the boundary displacements 
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A a = F(N - 1) [± h[ N - i) - El S li \~ im j ^ t W - V 



(3) 



(4) 



Let us introduce the dimensionless variables, Ki = ki/k, Si = Aj/Z, = Di/l and / = F/(klN), then Eq. 
the expression for the strain of spring i, viz. 



where 



e % = NfKi = Nf 



x? = {N-i)-s N (i-iy, 



X 



N 



E^ 3 «i(i-l)(^-i) 



= 1 fc 

We now express the total dissipated energy due to breaking of springs 



(5) 



(6) 



with properties X$ = —SnX^ 1 and Ki = 8(1 — Si/ri), where 0(x) = 1 for x > and 8(x) 



for x < 0. 



1 * 



i=l 

which can be rewritten in terms of the dimensionless variables in the form 

N 

kl 2 2 



(7) 



(8) 



i=l 



We further rewrite this equation in terms of the external applied force / using Eq. (0 , 



S{K} = - 



G{k} 



N 



(9) 



where G{k} is given by 



G{k} 



N 



3=1 



(10) 



Here, G can be interpreted as a global shear modulus of the lattice. It depends on the number of broken springs 
and the order in which the springs break. 

We define the order parameter for the system 



1 N 

-y 



(11) 



Then Eq. Q is expressed in the form 



F if ( P 

£ = - > Ki' 



2 tl IW) 



(12) 
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where G(0) is re-expressed as a function of the order parameter. As has been shown by Zapperi et al. the 
self-consistency requirement on yields the following integral equation (see Appendix) 

0=1- l 7 ^ tf(r)dr, (13) 
Jo 

where ip(r) is the probability distribution function for the system breakdown. 
The susceptibility for the system breakdown is then calculated according to 



^0 = #(y) 
df f '- y 2 il>{v)h<t,' 



y = ^ - (U) 



where 



/ . _ 1 G(0) + 0G(0) 

with G' = dG/d(j>. At the onset of breakdown, / = f c and = (f> c , the susceptibility diverges; this corresponds to 

n 20 e G(0 e ) 

yM = GWTW)' (16) 

with y c = / C /(^ C G(0 C )) 1 / 2 . We note that by specifying the form of the distribution function tp, Eq. iflfijl determines 
G(0 C ). In particular, for tp = c, where c is a constant, the differential equation llKit can be solved to give 

c 2 f 2 

G(0 C ) = 2' 

0c (00 - 0c) 

where 0o > C is an integration constant, which can be determined by computation. 
For f < f c , in the vicinity of f c , the susceptibility scales as 

x = w ~ (/c " /)_1/2 - (18) 

Similarly, the corresponding description for the order parameter is 

0-0c~ (/c-/) 1/2 . (19) 

Both relations Q18J1 and l|19ll are the mean-field theory predictions, see Appendix. 

The shear modulus G(0) is expected to obey, in general, a universal scaling law for ss 0, in the form 

G(0) = G O 6(0-1)0 T , (20) 

where Go is a material constant and r is a universal constant determined by experiment or computation. It is 
determined by re-arranging Eq. l|20"j) 

T=lim !2^Wl. (21) 

0^0 log |0| 

Finally in our analysis, we consider the concept of burst avalanche distribution [JJ. This concept has relevance to 
a system undergoing wear, for which the size of the avalanches can be used to estimate the size of the wear fragments. 
The size of a burst avalanche, s, in our system is defined as the number of springs that break simultaneously for a 
constant load The burst size distribution, Af(s), is generated by counting the number of bursts of size s that occur 
as the system breaks down. Hemmer and Hansen [jjj showed that for a large class of failure threshold distributions 
in the fibre bundle models, the bursts are distributed according to an asymptotic power law, 

lim^W*, (22) 



with a universal exponent, £ = 5/2. 
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FIG. 1: N parallel springs, connecting bodies I and II, 



subjected to a shear load F. 



III. COMPUTATION METHOD 



The mean-field theory presented above requires as input the shear modulus, G. Once the state of the system 
is known, i.e. the number and position of the ruptured links is known, the shear modulus of the system can be 
determined from Eq. Ijlflll . However, the position of the ruptured springs depends on the sequence in which they 
rupture. The sequence, in turn, depends on the local rupture thresholds and the applied force. Even for moderate 
system sizes (N ~ 1000) the number of possible sequences to failure become prohibitively large, hence a statistical 
method is needed to generate an average, or "most probable", shear modulus curve as a function of the order parameter 
<j). We take two different approaches to accomplish this; (a) a numerical scheme and (6) an analytic procedure. We 
compare the mean-field predictions with quasi-static simulations of the lattice breakdown. 



The probability that the spring i breaks first under a uniform strain field is p = 1/N, since the rupture limits 
are selected from the same distribution. Here the springs are subjected to a linear strain field, meaning that the 
probability that the spring i breaks first, can be written in the form 



where Wi is a weighting factor based on the local strains ef of the springs and the plus superscript indicates that only 
positive strains are counted, cf. Eq. 

The sequence of spring ruptures is then determined numerically by selecting a random number from a uniform 
distribution in the interval [0,1], and breaking the corresponding spring, see Fig- El The rupture probabilities, 
according to Eq. l12.'it . are then re-computed and a new rupturing spring is selected. The procedure is repeated until a 
complete path to total failure is established. The average shear modulus curve is then found by averaging over several 
such sequences. Finally, simple regression analysis is made to find a suitable polynomial or power law expression 
approximating the average shear modulus curve. 

A quasi-static method is used to perform the numerical simulations of the breakdown of the spring system (Fig.Q. 
The simulation is initiated by associating a rupture limit to each individual spring in the system. The limits are 
chosen from a predefined distribution. A full set of rupture limits is called a configuration. For each configuration, 
the computation starts with all links intact, Ki — 1, and the external load equal to zero, / = 0. Next, the load is 
incremented and the local strains are computed according to Eq. ©. The local strains are then compared to their 
respective rupture limits. If the local strain exceeds the rupture limit, then the link is broken, i.e. the local stiffness 
is set equal to zero. If a rupture has occurred, the local strains are re-computed and checked once more against the 
rupture thresholds, otherwise the computation proceeds with the next load increment. When only one spring remains, 
the load loop is interrupted and the computation proceeds to the next configuration. 



A. 



Numerical 



Pi = Wip, with u>i ~ 




+ ' 



(23) 
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Three different failure distribution functions are studied here; (i) a uniform distribution in the interval n £ [0;2], 

^( r ) \ otherwise 
(«) a Weibull distribution with shape factor v = 2 and location parameter A = 1, 

^(r) = |AU; (25) 
(Hi) a Gaussian (normal) distribution with mean value ■& = 1 and standard deviation a = 1, 



1 



2 /o„2 



tf(r) - — =e"^ . (26) 
crv27r 



The critical values, 4> c and / c , can be estimated directly from the mean-field analysis from Eq. iflfijl . For the uniform 
distribution this equation is easily solved, leading to a functional relationship between f c and 4> c - For the Weibull and 
normal distributions, a Newton-Raphson scheme is employed to find the load corresponding to a given value of the 
order parameter, by finding the roots of Eq. llfill . The result is a line in the (<fi c , f c ) plane, which can be compared 
with the critical values from the quasi-static simulations. 



P. 



FIG. 2: Illustration of how the sequence of breaking of springs is determined from a random number generated according to a 
uniform distribution, see Eq. 12,'ill . 



B. Analytical 

The material constant Go, cf. Eq. <20ll . is determined exactly when all the springs are intact, i.e. at <fi = 1. For 
this case, Eq. ifTfil) reduces to the following simple relation: 

1 + N 

GL i = G = . (27) 

Hence, as N — > oo, Go — > 1/4. The exact expression for G consists of terms like YliLi K ^ YliLi K i* an ^ 2i=i K »* 2 - 
With only one spring broken, denoted by j, the sums can be evaluated 

N 

= N-l, (28) 

t=i 

N 

= N 2 /2 + N/2 — j, (29) 

i=i 

JV 

^K.i 2 = N 3 /3 + N 2 /2 + N/6~ j 2 . (30) 



i=i 

Hence, the slope of G for cf) » 1 can be determined, once the position of the first spring rupture is known, viz. 



dG 



^ G(l) - G(l - 1/N) 
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We note that the first spring breaks according to a linear distribution, cf. Eq. 12 .'til . which is identified by j = 
int{2(JV — l)/9}. It corresponds to the average value of the linear distribution. Inserting this j into Eqs. ll29)l - (l30ll . 



then evaluating Eq. 1-Tlll and letting JV — > oo leads to ^ 
of G is found 



Consequently, the following linear approximation 



G +« = \+-k- (32) 

For the system under study, the JVth spring is always the last to break, since it is subject to compressive load. 
Moreover, the shear modulus is zero when only one spring remains intact. When only two intact springs remain in 
the system, denoting the next to last spring to break as m, the sums in G, cf. Eqs. p8 jl -l|30 fl . are evaluated: 

N 

t=i 

N 

^2k4 = N + m, (34) 



JV 
i=l 



i 2 = N 2 + m 2 . (35) 



Inserting Eqs. ll33|l - l(35|l in Eq. ifTf))) and using <j> = 2/JV, it is possible to carry out the limit analysis in Eq. 
once the position of penultimate spring rupture is identified. 

The springs closest to the lower side of the system (cf. Fig. ^| are subject to compressive stresses during the 
major part of the breakdown process. Assuming that the compressive stresses are sufficiently large, it is possible to 
identify the next to last rupture as m = JV" — 1. Carrying out the limit analysis leads to r = 3. Hence, the analytic 
approximation for G in the limit <\> ~ becomes 

G^o = \<t? ~ H 3 . (36) 
IV. RESULTS 

The results of our computations of the shear stress, the order parameter vs. the applied shear force and the burst 
size distribution are displayed in a number of figures. FiguresEb-fo show the shear modulus averaged over 100 rupture 
sequences of 2001 springs. The sequences were generated according to the procedure described in Sec. IIII Al This 
curve, henceforth, is called the average shear modulus curve. These figures also show the lower and and upper bounds 
of G(<p), respectively. The lower bound curve corresponds to the breaking of the highest strained spring and the upper 
curve to the breaking of the least strained spring. From these figures, it can be deduced that G(4>) depends linearly 
on (j) close to <j> = 1, while it has a scaling behaviour (power law dependence) close to = 0, cf. Eq. l|2"0j) . Note that 
the average and lower bound shear modulus curves coalesce as <j> — > 0, see Fig.|3>. This validates the assumption that 
spring JV — 1 is the second last to break and consequently that Eq. l(36|l is correct. 

Using regression, linear and power law expressions are fitted to the computed data. Here, the intervals S [0 + ; 0.25] 
and <p <E [0.75; 1.0] are chosen to determine the power law and linear law, respectively. Our analysis leads to the 
following expressions for G(<p), 

G^ ~ O.2170 3 1297 , (37) 
G ~i 0.22259 + 0.471540. (38) 

Figures show the analytical (Eqs. l-Kill and 1320 ) and numerical (Eqs. I|37ll - lj38ll ) scaling and linear laws, 

respectively, superimposed on the average shear modulus curve. As can be expected, the approximated curves deviate 
for <f> in the mid-range (<J> ~ 0.5). 

A system size of JV = 10001 springs and a load increment of A/ = 0.01/JV are used in the numerical simulations 
(the quasi-static method discussed in the foregoing section) presented here. To gather enough data for statistical 
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predictions, each simulation is repeated for M = 100 different configurations of rupture thresholds r. Values from the 
last load step before global breakdown are taken as critical values. Figure displays the simulated critical values 
of the order parameter, <j> c and the shear force f c , for the uniform, Weibull and normal distribution functions defined 
in Eqs. I|24ll - lj26|) , The three distributions display similar results with respect to scatter in f c and cj> c , see Figs.|3>-c. 
Using these results (average values of G((p c ), <f> c and / c ), we have determined <f>o in Eq. I|17J1 for the uniform failure 
distribution given by Eq. We found (j>o = 1.042. 

Next, we focus on the relationship between the scaling variables <j> and /. Figures Eb- c show the computed order 
parameter, </>, plotted against the applied load, /, for the three distributions using the quasi-static method. The 
data presented in the figures are averaged over 100 configurations of rupture threshold. The displayed {cf>, /)-plane 
is divided into rectangles. The arithmetic mean of the values confined within each rectangle is used to represent the 
data in the figures. The number of points within each rectangle is utilised as weighting factors to fit a scaling law 
expression Eq. Ill9fl to the simulated data. The scaling form of the order parameter predicted by the mean-field theory 
Eq. Q19J1 is also depicted in Figs. UJi-c. The simulated data support the validity of the mean-field analysis. 

Finally, the important concept of burst distribution is explored. For a system undergoing wear, the size of the 
avalanches can be used to estimate the size of the wear fragments Q. The size of an avalanche or burst, s, is defined 
as the number of springs that break simultaneously for a constant load • The burst size distribution is generated 
by counting the number of bursts of size, N(s), occurring as the system breaks down. Figures 03- c show the burst 
distributions for the three chosen rupture threshold distributions for the damage model considered here. The three 
figures display similar behaviour. However, our simulations indicate that the Weibull and normal distribution give a 
higher amount of large bursts (s > 10). The figures also include a line with a slope —5/2 (log-log scales), which is 
the theoretical result found for a fibre bundle model with global load sharing pj], ■ The results indicate that this 
scaling law has relevance for the system studied here. 




(a) (b) 

FIG. 3: Numerical simulation of G(<j>) averaged over 100 randomly ordered series of ruptures of 2001 springs, solid line. The 
dotted and dashed lines depict the lower and upper bounds of G, respectively, (a) and (6) show linear and logarithmic plots, 
respectively. 



V. DISCUSSION 

The results presented in the preceding section show that the average shear modulus is linearly proportional to the 
damage in the beginning of the breakdown process and has power law characteristics at end of the process, see Fig.EJ 
Consequently, the shear modulus can be approximated with simple relations, which do not require knowledge of the 
path to failure, i.e. the order in which the springs break. The critical force and critical order parameter (damage) 
can be estimated analytically using mean-field theory once the shear modulus has been determined, provided that the 
distribution of the disorder of the material is known. Here, it is assumed that the disorder is quenched in the system 
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FIG. 4: Approximate computations of the shear modulus G(4>). The solid and dotted lines display numerical and analytical 
approximations, respectively. The broken line shows the curve found from averaging over several sequences of spring ruptures, 
see Fig. El 



and does not change during the breakdown process [lfij . 

In the preceding sections, we showed that for 4> w 0, the shear modulus G scales as G ~ (j) T , with r = 3, which we 
obtained by analytical calculation, while numerical simulations led to t w 3.13 in the interval S [0 + ;0.25]. These 
results may be compared with properties of other condensed matter systems, e.g. the macroscopic conductance in a 
resistor network or the elastic modulus of gels in polymer-system j2jj]. For lattice dimensions d > 6, which corresponds 
to a mean-field approximation, de Gennes found that r = 3. This result also conforms with a recent calculation by 
Xing et at of the scaling behaviour of shear modulus in polymer-system j^Jl . Their renormalisation group analysis 
shows that r = 3 — ^-e, with e = 6 — d, which is identical to the conductivity exponent of the random resistor network 
calculated by Harris and Lubensky [22j. 

The constitutive relation for damage, <f> vs. /, for all the three statistical distributions (Figs. HJj-c) indicates that <j> 
decreases slowly with the increase in /, then a sharp (discontinuous) drop in </> occurs, when the shear force reaches a 
threshold value f c . This sudden transition in the order parameter and the fact that the scaling form given by Eq. I|19JI 
describes the approach to f c so well (/ < f c ), indicates that material breakdown is a first-order phase transition, e.g. 
similar to a spinodal decomposition, see Appendix. 

The burst distributions for the system using the three statistical distributions are displayed in Figs.[J|i-c. Here the 
size distribution of avalanches is integrated over the load. For a fibre bundle model (FBM) with global load sharing 
(the load is evenly distributed over the remaining fibres) the expected number of avalanches A/(s), where s is the 
size of the avalanche, scales as N ~ s _? , with £ = 5/2 Although this scaling law holds for a large class of 

threshold distributions, this is not always the case j23|. For small avalanches, s < 10, a regression analysis, using 
the numerical data shown in the figures, gives an exponent close to the FBM value, 2.5 ± 0.2. It may be worth to 
compare our results with that for the random-field Ising model (RFIM) predictions, which intends to describe the 
behaviour of many magnetic materials subjected to external magnetic field H. For this system, the magnetisation 
changes through nucleation and motion of domain walls. The motion is discrete, i.e. there will be magnetisation 
bursts corresponding to reorganisation (or avalanche) of a domain of spins which span several decades of size [24l |. 
There will be plenty of small avalanches of spins and fewer and fewer avalanches of larger and larger sizes . The 
power law dependence for the probability of having an avalanche of a given size has been determined [24[ ■ It is shown 
that the size distribution of all avalanches that occur in one branch of hysteresis loop (for H from — oo to +oo) scales 
in the form: Mint ~ s~ e J 7 (s, r) , where T is a scaling function and r is a "disorder" parameter j2^|. The mean-field 
computation of Perkovic et al. 26] gave 9 = |, which is close to the corresponding value predicted by FBM. 

The breakdown of the system that we considered is similar to that of FBM. The process of breakdown is envisaged to 
proceed by the formation of micro-cracks. The micro-cracks are nucleated according to a linear distribution function, 
which is proportional to the displacement field. The cracks grow up to a critical size and coalesce to form the final 
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FIG. 5: (a) The critical scaling variables 4> c versus f c for 100 configurations, computed by the quasi-static method, of uniformly 
r\ € [0;2], Eq. I|24fl . circles, Weibull (y, A)=(2,l), Eq. 125(1 . squares and normal (#, cr)=(l,l), Eq. H2tifl . diamonds, distributed 
rupture thresholds of N = 10001 springs. (6) and (c) display normalised histogram distributions of C and f c , respectively. 
Solid, dotted and dashed lines refer to the uniform, Weibull and normal distribution, respectively. The abscissae have been 
rescaled with the mean value, and the standard deviation, a, to enable comparison between results for the three rupture 
threshold distributions. 

crack. We found no indication that a single crack takes over and dominates the breakdown process. This behaviour 
was also found in a two dimensional mesh subject to a uniform stress field 



VI. CONCLUSION 



We modelled an interface layer by N parallel elastic springs that connected two rigid blocks. We loaded the system 
by a shear force / acting on the top side. The springs were assumed to have equal stiffness but were ruptured randomly 
when the load reached a critical value f c . The state of the material damage for the system, was characterised by 
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FIG. 6: The damage order parameter <j> versus the shear force / averaged over 100 configurations, computed by the quasi-static 
method of (a) uniform n 6 [0;2], (b) Weibull (v, A) = (2,l), and (c) normal cr) = (l,l), distributions of rupture thresholds for 
N = 10001 springs, respectively. The lines represent the mean-field scaling behaviour. 



an order parameter, <fi, varying between zero (complete breakdown) and one (intact). The shear modulus for the 
system was numerically simulated over the entire range of <f> and it was found that G^-o ~ 313 and G^-i <~ <f>. The 
linear behaviour of G((j>) close to <fi = 1 was confirmed analytically. Moreover, our analytical calculation showed that 
G^-o ~ <j) 3 - We calculated <fi vs. / numerically and determined their values at the onset of breakdown, i.e. around 
(4>c,fc), and also found their variation depending on the failure distribution function selected. We noted that there 
was a discontinuous drop in 4> around / = f c . The mean-field theory predicted this behaviour for / ~ / c , leading 
to the conclusion that breakdown behaves like a first-order phase transition. Finally, we calculated the burst size 
distribution during rupture and found that the system behaves according to the predictions of the fibre bundle model. 

The model presented here can be extended to treat a more realistic representation of wear induced loads by modelling 
a two-dimensional lattice comprising a network of blocks connected to each other and to a rigid indenter by elastic 
springs. 
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(c) 

FIG. 7: Burst size distributions for 100 configurations using (a) uniform rj £ [0;2], (6) Weibull (v, A) = (2,l) and (c) normal 
<r) = (l,l), probability distributions of rupture thresholds for N = 10001 springs, respectively. The straight lines have a slope 
of -2.5, see Eq. J22j. 
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APPENDIX: MAGNETIC SYSTEM AND LATTICE DAMAGE: AN ANALOGY 

In magnetic language (Ising model) , let us envisage that the lattice sites confine "spins" which may have only values 
(T.j = +1 and (Ji = —1. These spins interact with the nearest-neighbours on the lattice with a constant interaction 
energy Jy, see e.g. In particular, we consider the random field Ising model (RFIM), where at each site i there 
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is a random field gi and there is an external field H acting on the system, the Hamiltonian is given by [2£ 



N N 



<i3> 



The system transforms from negative to positive magnetisation as the external field H is swept upward. Supposing 
that each spin a t always gets aligned with the effective local field and Jy = J, the Hamiltonian IIA.lll is simplified to 



N N 



Wj = Jo,- +# + <?;). (A.2) 



In the mean-field theory of Weiss (0J V 29]), one assumes that every spin is coupled to all N — 1 other spins with 
a coupling constant J/(N — 1). The effective field acting on site i is now expressed as hi = Jzj\4/(N — 1) + H + gi 
where z is the lattice coordination number and M. is the magnetisation, the order parameter, defined as 

N 

M = °'j = N+ — N- , (A.3) 

3 

where N + = (N + j\4)/2 and 7V_ = (N — j\4)/2 are the number of up-spins and down-spins, respectively. Moreover, 
spins with gi < —(Jzj\4/(N— 1) + H) will point downward, others will point upward. Since gi is a stochastic variable, 
the self-consistency j2^| provides 

/— Jzm — H 
P(g)dg, (A.4) 
-oo 

where m = M./N is the magnetisation per site (N » 1) and p{g) is the probability distribution for the random field 
g. The magnetic susceptibility is calculated by differentiating Eq. l!A.4|l with respect to H, viz. 

X= dH = l- 2jl{x) ' (A - 5) 

where x = — Jm — H and we put z = 1 for convenience. Note that at x c = p _1 (0.5/J) the susceptibility diverges. 
Taylor expansion around x c yields, 

X - J^ 1 + 0.5J- 2 ^ 

p'(x c )(x - x c ) + hp"(x c )(x - x c ) 2 + ■ ■ ■ 



Near x c — —Jm c — H c , we have 



p'[X c )(X - Xc) 



Re-writing Eq. l|A.7Jl in terms of m and H, 

dm A 



dH J(m - m c ) + (H - H c ) ' 



(A.f 



with A = 0.5J 2 / p'(x c ). Integrating Eq. flA.Sjt . we find the equation of the state in transcendental form for H < H c 
and, close to H c , viz. 

m-m c = - J Blog[l + ^ ^-2 ], (A.9) 

with B — J A. Again series expansion of Eq. IIA.9J) around m c and H c results in the dominating contribution, 

m-m c ~ (Hc-H) 1 / 2 . (A.10) 
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The susceptibility in this region (H ss H c ) scales as 

X-iHc-H)- 1 ' 2 . (A.ll) 

The scaling behaviour for the magnetisation, Eq. ijA.lOl! and the susceptibility Eq. (IA.11I) . close to H c for H < H c , 
resemble that of the spinodal decomposition in the mean-field theory of nucleation (lf| . 

For the lattice damage model, upon analogy with RFIM, we write the following Hamiltonian 

N N 

U D = ^2 eKiK i + ^ W '(M ~ Pi)> ( A -12) 

<ij> i 

where e is the inter-spring interaction parameter, fj, is the effective chemical potential and pi is a random field 
variable. To make contact between the lattice damage and magnetic system, we write: <x; = 2ni — 1. The lattice 
damage Hamiltonian, Eq. (|A. 121 . expressed in terms of spin variable Ci, becomes 



^ JV N N 



4 

<ij> %— 1 i— 1 

where we have left out the constant spin-independent terms. Comparing Eqs. i|A,13Jl and l|A.lj) . we can relate the 
variables in the magnetic model and the lattice damage model as follows 

J=-e/4; H = -(e + n)/2; 9l = Pl /2. (A.14) 
Let us first rewrite Eq. l|A.12j) in the form 

N N 

Hp = Kj (ey~] «j + M - Py i (A. 15) 

i 3 

then introduce the order parameter $ for the lattice damage model 

N 



$ = ^Ki = N -q. (A.16) 



Here, (N — q) and 9 are the number of intact and broken springs, respectively. In the mean-field theory Eq. ljA.1 
is replaced by 



N T 



ftD=£>i(^7 + /i-Pi), (A-17) 
Fashioned after Zapperi et al. 0, E| , we rewrite Eq. JjjJ in the form 



„2 



2 fe V*G($) 



(A.18) 



where G($) is re-expressed as a function of the order parameter $. Now comparing Eqs. l!A.17j) and we 
identify 



f 2 



where we placed 4> — ®/N. 
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The relation between (f> and the magnetic order parameter is m = 20 — 1 . Using this relation and Eq. l|A.14j) , we 
can write Eq. l|A.4jl as 



p(g)dg = 1 - - / 

-oo L JO 



2e<j>+n 



p(p)dp, 



(A.20) 



where we placed z — 2 and considered only pi > 0. By invoking Eq. l|A.19|l . we find 



cf> = 1 j * W p(r)rdr. 

2 Jo 



(A.21) 



Similar calculations as made on the magnetic order parameter m yield the scaling relations for the lattice damage 
order parameter 4> and the susceptibility x, cf. Eqs. ijA.lflll and l|A.llj) : i.e. close to 4> c for / < / c , we have 



~ (fc-f) 1/2 , 

*=f -(/-/)-". 



(A.22) 
(A.23) 



where / c is the critical force for the lattice breakdown. Both relations l|A.22l) and (IA.23|) are the mean-field theory 
predictions. 
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